NUMERICAL SIMULATION OF GLUEY PARTICLES 
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Abstract. We propose here a model and a numerical scheme to compute the motion of rigid particles 
interacting through the lubrication force. In the case of a particle approaching a plane, we propose 
an algorithm and prove its convergence towards the solutions to the gluey particle model described 
in We propose a multi-particle version of this gluey model which is based on the projection of 

the velocities onto a set of admissible velocities. Then, we describe a multi-particle algorithm for the 
simulation of such systems and present numerical results. 

Resume. Nous proposons ici un modele ainsi qu'un schema numerique afin de resoudre le mouvement 
de particules rigides en interaction a travers la force de lubrification. Dans le cas d'une particule a 
I'approche d'un plan, nous proposons un algorithme et montrons sa convergence vers le modele de 
particules visqueuses decrit dans [15]. Nous proposons une version multi-particules de ce modele qui 
est basee sur la projection des vitesses sur un espace de vitesses admissibles. Ensuite, nous decrivons un 
algorithme multi-particules pour la simulation de tels systemes et presentons des resultats numeriques. 



1. Introduction 

Slurries, lava's flows or red cells in blood are systems made of rigid particles embedded in viscous 
fluids (if we consider as a first approximation that red blood cells are rigid). Such systems can also 
be found in industry: concrete, paper pulp or some food industry products. These systems present 
varieties of noticeable rheological behaviours, whose study has been the subject of a great amount 
of researches, with contributions coming from engineering, chemistry, physics or mathematics. The 
basic problem is to predict macroscopic transport properties of these suspensions - viscosity, settling 
velocity - from microstructures, that is to say, from the interactions between particles and from their 
spatial configuration. 

In case of dilute suspensions, theoretical results come from neglecting near field interactions. For 
example, in 1906, Einstein proposed an asymptotic formula for the apparent viscosity of dilute sus- 
pensions [S]. In that case, apparent viscosity only depends on the solid volume fraction. Unfortunately, 
agreement between such asymptotic results and experiments generally fails as soon as the solid fraction 
reaches a few percent. For higher solid fractions, near field interactions can not be neglected anymore 
and it becomes essential to take them into account. Note that, studying the behaviour of neighbouring 
particles is of great interest, not only to understand the behaviour of dense suspensions, but also to 
study the fluid/particle system of equations modelling suspensions of particles. Indeed, existence of 
solutions to these equations has been proved as long as the distance remains stricly positive (see for 
example [H EH ES] ) . Global weak solutions have also been constructed in [6l[22], supposing that solids 
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stick after contact. However, nothing is said concerning the possibility that such a contact may occur 
in finite time. A good understanding of near field interactions is therefore necessary to study more 
precisely these systems. 



These interactions between solids embedded in a viscous fluid are due to lubrication forces: for the 
solids to get very close, the fluid must be evacuated from the narrow gap between them, which creates 
a force penalizing their relative motion. This force is singular in the distance and this singularity is 
sufficient to avoid contacts. Indeed, it has been proved in [9] that in two dimensions, a smooth particle 
embedded in a viscous fluid following Navier-Stokes equations can not touch a plane in finite time. 
This behaviour can be recovered from the asymptotic expansion of the lubrication force, available for 
a Stokes fluid in three dimensions (see [T] for example): 

(1) Flub ~ -QTTfir^^, 

q 

where fi is the viscosity of the fluid, r the radius of the particle and q the distance between the particle 
and the plane. Indeed, using this first order approximation, we can write the Fundamental Principle 
of Dynamics for a particle of mass m submitted to an external force /: 

(2) mq{t) = —GiTfir^- + mf{t), 

q 

and the fact that the maximal solution to this ODE is global and never goes to zero (contact) in finite 
time comes from the Cauchy-Lipschitz theorem. Similarly, in case of a fixed sphere of radius ri and 
another sphere of radius r2 moving at velocity V along the axe of the centers, the first term of the 
developpement of the lubrication force exerted on the moving particle is (see [2]): 

(3) Flub ~ -67r/i 
and no contact can occur in finite time. 



This force, while acting at microscopic level, can be very important for the macroscopic behaviour 
of the global system, especially in case of high density of particles. Even for Stokes fiows, it induces 
complexity and nonlinearity. This complex link between microscopic and macrosopic levels makes 
it difficult to obtain theoretical results and studying these systems requires numerical simulations. 
In order to obtain relevant simulations for dense suspensions, the lubrication force has to be taken 
into account with accuracy. However, direct numerical simulation induces space discretization which 
makes it difficult to solve accurately the fiuid in the narrow gap between neighbouring particles. As a 
consequence, numerical contacts can be observed in such simulations and physical reasons as well as 
numerical robustness make it necessary to develop specific technics to deal with these contacts. 



A first idea to solve this problem is to search for a strategy allowing an accurate computation of 
the lubrication forces. In [TO], a method based on local refinements of the space and time meshes is 
proposed, so that the lubrication force in the interparticle gap is taken into account with accuracy 
and prevents overlappings. However, the number of refinements needed is not known a priori and the 
method can become computationally heavy. Consequently, less time-consuming methods have been 
developped. Some of them consist in adding a short range repulsive force (see [3 [21] or [27]). In [17] 
a minimizing algorithm is used to impose a minimal distance between the particles, while in [23], 
the particles are allowed to undergo slight overlappings and an elastic repulsive force is added when 
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such overlappings are detected. All of these methods ensure numerical robustness but introduce new 
parameters and do not take into account the underlying physics. Another approach is to use inelastic 
collisions. This idea has been proposed in [11] in order to impose a minimal distance between the 
particles. In [18], a scheme for inelastic collisions, based on a global projection step of the velocities, 
has been developped for granular flows and makes it possible to handle lots of particles. This scheme 
has been coupled with a fluid/particle solver in [13], to avoid contacts. More physical strategies, taking 
the lubrication force into account, have finally been proposed. Each of them relies on the asymptotic 
developpement of the lubrication force ([3]). In [3l [20], it is shown that these lubrication forces are 
solution to a linear system. They are computed at each time step and added to the simulations. 
Unfortunately, this leads to stiff systems and, whereas it better takes into account the underlying 
physics, contact problems still occur because of the time discretization. In [16j a method is proposed 
to stabilize this problem by computing accurately sensible quantities such as the interparticle distances. 
However, a projection step is still needed for big time steps, in order to avoid overlappings. 

The purpose of this article is to propose a strategy dealing simultaneously with contacts and lubrication 
forces. We restrict ourself here to the study of a gluey contact model without taking the surrounding 
fluid into account. This model is based on the gluey particle one described in [19]. We propose an 
algorithm for this particle/plane model and prove its convergence. Then, we generalize it to the multi- 
particle case. The numerical strategy is to combine the algorithm given for the plane/particle case 
with the scheme proposed in [18j for granular flows. While programming this multi-particle algorithm, 
we watched out for dealing with contacts efficiently in order to manage to simulate collections of many 
particles. Numerical simulations for few thousands of gluey particles are presented in the last section. 
An example of coupling with a fluid/particle solver is given in section [2741 in the particle/plane case. 



2.1. The gluey particle model. We consider a three-dimensional spherical particle moving perpen- 
dicularly to a plane (See Fig. [T|). Its velocity and radius are denoted by V and r respectively. Its 
distance to the plane is q. 



2. Single particle above a plane 




Figure 1. Notations. 
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The gluey particle model has been proposed in [19]. It describes, from a macroscopic point of view, 
the behaviour of the system near contact. It is built as the vanishing viscosity limit of the lubrication 
model ([1]) and relies on two states, glued {q = 0) or unglued {q > 0). These states are described by 
a new variable 7 which stands for an adhesion potential: the more 7 is negative, the more the solids 
are glued. 

We denote by / =]0, T[ the time interval. The unknowns q and 7 belong to the following functional 
spaces: 

q e , q G BV{I) , 7 G BV{I), 

and the initial conditions are: 

g(0) = g°>0, g(0)=u°, 7(0) = 0. 

In order to be able to generalize the model to the multi-particle case, we use the following second 
order ODE formulation given in |19j : 



(4) q{t+) = nc„,(i)(zr ), 

(5) mq = mf + X in M{I) = {Ceil))' , 

(6) supp(A)c{t, q{t)=0}, 

(7) 7 = -A, 

(8) (?>0, 7<0, 



where Cq^^{t) is the set of admissible velocities at time t: 

{0} if ^{t-) < 0, 
R+ if^(t-)=o, qit) = 0, 
R else . 

Remark 2.1. In this formulation, q and 7 are supposed to be in BV{I). In order to alleviate the 
notations, their differential measures have been denoted by q and 7 respectively. 

The behaviour of the solutions to this problem is the following. By ([5]) and ([6]), g is solution to q = f 
while there is no contact (q > 0). Suppose a collision occurs at time to, we have (/(tg ) < and 7(to ) = 
0. Then Cq^y{to) is and ([4]) gives ^(i^) = 0. By ([5]), we obtain that, in the sense of distributions, 
A identifies to the dirac mass at time to weighted by the velocity jump — g(to )) = ~iT^Q{to)- 

This, together with ([7]) finally gives that 7 is initialized to the value mqit^) < 0. From then, while 
7 remains strictly negative, Cg^^ is reduced to {0} and, combining this with (jj]) gives that there is 
adhesion between the solids (q = 0). During this adhesion, q is zero and therefore, ([7]) associated 
to ([5]) gives 7 = mf. By definition of Cg^-y, the particle is allowed to take off when 7 is back to zero. 
An example of such a behaviour is given in figure [2] 



Cq,y{t) — 
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Figure 2. Example of solution to the gluey particle model. 

Remark 2.2. The additional constraints ([8]) are necessary. Indeed, suppose that t\ is an unsticky 
contact time (^'(ti) = 0, 7(tj~) = 0, = 0). If the force is negative after this instant and if we do 

not impose (7 > 0, then 7 = and q{t) = J^^ f{s)ds is a solution to the problem and the particle can 
enter the wall. Similarly, if the force is positive and if the constraint 7 > is not imposed, q = and 
7(t) = ?Ti Jj* f{s)ds is a solution and 7 can become strictly positive. 

Before proposing an algorithm to compute the solutions to this model, we make a few remarks about 
its interpretation and its physical relevance. 

Remark 2.3 (Physical interpretation). As already mentioned, a smooth particle embedded in a 
newtonian fluid never touches the plane in finite time. In the context of the gluey particle model, the 
variable q can be seen as a macroscopic distance between the solids: it is equal to zero as soon as the 
solids are near contact. The new variable 7, which is obtained as the limit of 7^ = (j.ln{q), stands 
for the microscopic distance. To understand the behaviour of the gluey particle system presented on 
figure m one can consider a rigid ball falling on a table coatted with a viscous fluid like honey. When 
the particle reaches the layer of fluid, it instantaneously sinks in it and the depth it reaches is linked 
to the impact velocity. From then, the ball is glued to the layer of fluid, the macroscopic contact 
begins, q is set to zero and 7 stores the impact velocity. As long as it is pushed, the particle sinks 
deeper in the fluid and gets closer to the plane (7 decreases). Then the particle is pulled. From that 
moment on, it smoothly moves back from the fluid (7 increases) and comes unstuck from the layer 
of fluid when the pulling forces have balanced the impact velocity and the pushing forces (7 reaches 
zero). Note that, from ([5]), A can be interpreted as an additional force, exerted by the plane on the 
particle, in order to satisfy the constraint It follows from ([6]) that the plane is allowed to act on 
the particle through this force only if they are in macroscopic contact. 

Remark 2.4 (Radius). This gluey particle model is built in [19] as the vanishing viscosity limit of 
the lubrication model ([T]) where each constant except the viscosity is taken equal to 1. Taking all 
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constants into account leads to define 7 as the limit of 7^ 
evolution becomes 



^■n^ln{q) and the equation governing its 



(9) 



7 



The larger r is, the less the microscopic distance 7 varies (the more it is difficult for the particle to 
move). Note that, provided we are only interested in the macroscopic trajectory q of the particle, the 
previous model (Ill)-(l8|) was valid for any radius: these trajectories only depend on the sign of 7 (and 
not its value) which is independent of r from ([9]). 

Remark 2.5 (Viscous or not viscous ?). Since this model is built by letting the viscosity go to zero, 
one may wonder whether it models viscous fluids or not. To answer this question, we consider the 
same experiment as in figure [2] (pushing untill time 2 and then pulling) for a particle falling on a 
plane coatted with different viscous fluids. On flgureO we compare the trajectory given by the gluey 
particle model to the trajectories computed for these systems where the viscous fluid layer is modeled 
by (121) • Of course, trajectories converge to the limit model when the viscosity goes to zero. We also 
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Figure 3. Comparison gluey particle model / layer of viscous fluid. 



observe that, from a macroscopic point of view, as long as we are interested in hitting and unsticking 
times, the limit model seems to agree with all trajectories. As a matter of fact, from this point of view, 
what is important is not the viscosity of the fluid but whether it is viscous or not. However, the main 
difference between the trajectories is the minimal distance reached by the particle. Actually, small 
viscosities induce very small distances and, consequently, the system can reach a domain wherein the 
initial lubrication model ([2]) is no longer valid. To conclude, the gluey particle model shall be employed 
to represent the macroscopic behaviour of very viscous systems for which distances are not too small. 
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2.2. Numerical algorithm. We propose here an algorithm for problem dH)-®. Let h = T/N be 
the time step. The problem is initialized to > 0, G R and 7*^ = A*^ = 0. We denote by q^, u", 7" 
and A"" the computed values of q, u, 7 and A at time t". We define /" by /" = -j^ f{s)ds. We 
have to compute g"""*"^, u""*"^, 7""*"^ and A^^^. 

In order to compute li""*"^ and A"'*'^, we define the discrete counterpart of Cq^y{t^) the following way: 

Kiq"", 7^^) = {v , + /it; > 0} if 7*^ = 0, 
7") = {v , + /it; = 0} if 7" < 0. 
K{q^, 7") is called the set of admissible velocities at time t". The collision law ^ and the Fundamental 



Principle of Dynamics ([ 

„«+l/2 



,n+l 



then become, 



n+l 



,n+l/2 



mm — 

«G-fsr(ij",7") 2 



V — u 



n+1/2 



where {v,w)m, = {mv,w). Note that the velocity the particle would have at time t"^"*^ if 

there were no plane, -u""*"^ is the projection of this a priori velocity on the set of admissible velocities 
K(q"^, 7") for an adapted scalar product. From this projection step, arises a Lagrange multiplier, 
denoted by A""*"^ (positive if 7" > 0), and such that 



71+1 



u 



This can be rewritten as 
(10) 

which is a discretization of ([5]). 



u 



771- 



mr + A 



n+l 



Then, 7"+^ is given by an explicit Euler discretization of ([7]), 



This equation is valid while 7"+^ is negative. If it becomes strictly positive, it means that the particle 
has taken off at a time t* In that instance, 7""'"^/m has integrated the force on 

instead of u^~^^ which was fixed to zero. Therefore, in that case, we modify u""''^ and 7"+^ the following 
way : 



if 7"+^ > , u 



n+l 



7" /m and 7 



n+l 



0. 



Finally the position is given by 



f +1 = + /IU"+^ 
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To sum up, the algorithm is the following : 

Algorithm 2.6 (Particle/plane). For all n > 0, let , vP , 7" and A" he given. We define 

,Ti 4-1 

h Jtr> 



f{s)ds. 



[1) Computation of the a priori velocity, without taking the lubrication force into account 



u 



n+l/2 _ n 



(2) Projection of the a priori velocity onto the set of admissible velocities, 



2 1 
= min — 



V — U 



n+l/2 



,n+l 



where K{q,j) = {v , q + hv > 0} if ^ = 0, 
K{q, ■j) = {v , q + hv = 0}ifj<0. 
From this projection step, we obtain A"^"'^. 

(3) Updating ofj, 

(4) Modification if unsticking, 

if 7"+^ < 0, n"+i = and 7"+^ = 7"+\ 

if 7""^^ > 0, u"+i = 7"+Vm' o.nd 7"+^ = 0. 

(5) Updating of q, 

Q" 



Remark 2.7 (Coupling with fluid simulations). This algorithm simulates collections of gluey particles. 
Let us now suppose that the particles are embedded in a viscous fluid. To make simulations taking 
the lubrication force into account, a splitting method can be used to couple a fluid/particle solver 
with the gluey particle algorithm. We denote by and p" the velocity and pressure fields into the 
fluid at time t"". Let S be any fluid/particle solver: from u", q^ and S computes the a priori 
velocities of the particles, without taking the lubrication force into account carefully. To couple the 
two algorithms we propose to modify step ([T|) of algorithm 12.61 writing: 

,n+l/2 



U 



2.3. Convergence result. In this section, we establish a convergence result for the proposed algo- 
rithm. To begin, we rewrite problem ([I])-® as 

mq + ^ = m ^^(0) + j f{s)ds^ , 

(11) { q>0, 7<0, 97 = 0, 

q{0)=q^>0, g(0)=7X°, 

which is formally equivalent to the previous one (see [19] ) . 

We recall that h is the constant time step. We denote by qh the piecewise affine function with 
qh{t^) = q^- Similarly, 7/1 is the piecewise affine function with 7/1 (i") = 7"- We denote by Uh the 
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derivative of qh, piecewise constant equal to u"^^ on Finally, we define = — piecewise 

constant. Note that, due to step (jl]), Xh is generally not equal to A"+^ on We will denote 

by A"^^ = — — 7")//i its value on this interval. If the particle does not take off between times 
and t""*"^, no modification is made during step ([3]) and we obtain A"'"''^ = A"^^. The convergence 
theorem is the following: 

Theorem 2.8. Let f be integrable on I =]0, T[. When h goes to zero, there exists subsequences, still 
denoted by {qh)h, {uh)h, {\)h and {■^h)h, Q G nC{I) and 7 G BV{I) such that 

Uh — > u in L}{I), 

qh — > q in W'^'^{I) and L°°{I) with q = u, 

Xh — ^ A in M{I), 

7/i — > 7 in L^{I) with 7 = —A, 

where (^,7) is solution to f77]) . 

Remark 2.9. Non-uniqueness for the limit problem (see [19] for counter-example) prevents from using 
the standard approach based on consistance and stability. Consequently, we use compactness methods 
and obtain convergence up to subsequences. However, in case q has a finite number of zeros, the limit 
model admits a unique solution and therefore the convergence of the algorithm to problem (jlip is 
proved. Moreover, in that case, it can be shown that dl])-® and (llip are equivalent, in the sense that 
a solution to one of the problem is also solution to the other (the demonstration of this result can 
be found in [14J). Consequently, under the a priori hypothesis that q has a finite number of zeros, 
theorem 12.81 shows that algorithm 12.61 converges to Q-dHl). For example, this hypothesis is verified if 
the external force / changes of sign a finite number of times. 

Proof of theorem 12.81 

To begin, note that a discrete form of the Fundamental Principle of Dynamics ([5]) is verified: 

ra+l _ n 

(12) Vn, m- -^ = m/" + A"+^ 

h 

Indeed, in case the particle does not take off between times t" and t^^^, the equality follows from (jlOp . 
together with A"^^ = A"^^. If the particle takes off, it comes from (jlOp and step H] of the algorithm. 

The proof will be devided into 4 steps. 

(1) Convergence of qh and Uh 

Lemma 2.10. {uh)h is bounded in L°°(I). 

Proof : In case the particle does not take off, the projection step ([2]) gives 

|^n+l| ^ |-n+l| < |^n+l/2| < |^n| ^ ^|^n|_ 

In the other case, it can be proved that m""'"^A""'"^ < and combining this with (jl2p gives the 
same result. By summing up all these inequalities we obtain 

K+i|<|nO|+ Tl/I, 
Jo 
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and the result follows from definition of u^. □ 



Lemma 2.11. {uh)h is bounded in BV{I). 

Proof : By lemma [27101 the result will follow provided we prove Var(n/i) is bounded indepen- 
dently from h, where 

N-l 

VarK) = -n"|. 

n=l 

To check this, we first split the sum and consider the sums between indexes pi and ni where 
t^^ and t"^ are successive unsticking times (See Fig. H]): 

ni — 1 
n=pi 

The total variation of is made of a sum of such terms. 

• • • 

* • • — I • • • I — > 

f" iPi t"" t"i 



Figure 4. Proof of lemma [2TTT] : notations. 

The idea behind the above decomposition is that, at each unsticking time t^^ , the velocity of 
the particle is small and that its variations over [t^^ , [ only depend on the integral of / over 
the same interval. These terms will be summed up to obtain a bound on the total variation. 

More precisely, the bound for Var[jpi jni[(n/i) can be found by analysing each jump — 
li""!, paying attention to what happens at time (hitting time, sticking time, unsticking 
time). For the sake of readibility, details of the computation are skipped here, they can be 
found in [14]. We find 

Yar[tPi ^t'n[{uh) < 4 / \fis)\ds. 
Summing up all these contributions and the bounding terms, we obtain 

Yariuh) < u'^ + 8 [ \fis)\ds, 
Jo 

and Var(ti/i) is bounded independently from h as required. □ 

Lemma [2.111 together with the compact embedding of BV{I) in L^{I) gives (up to a sub- 
sequence) 

Uh — > u m L^{I) with u e BV{I), 
Qh — > q in W^'^{I) with q = u. 

Uniform convergence of qh to q then follows from the continuous embedding of W^'^{I) in 
L~(/): 

qt^qm L^{I). 
Finally, since q^ is positive, we have q>Q everywhere. 
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(2) Convergence of 
Lemma 2.12. {\h)h is bounded in L^{I)- 
Proof : By ([l2]) and the fact that = we get 

T 

\Xh\ < mNai{uh) + 

The result follows by combining this with lemma [2.111 □ 

By lemma [2.12| {\h)h is bounded in A^(/), which implies that there exists a subsequence 
and A G such that 

A,, ^ A in M{I). 

Moreover, combining lemma [2. 121 with 7/j = —\h it comes that {'^h)h is bounded in BV{I). 
This, together with compact embedding of BV{I) in L^{I), implies that there exists a subse- 
quence and 7 G BV{I) such that 

(15) 7/i — > 7 in L^{I) and a.e. 

Since 'jh is negative, it follows from this convergence result that so is 7. Finally, since ■jh = —\h, 
we can check that 7 = — A in 

(3) Continuous FPD 

We are now going to prove that mq + 7 = m ^q'(O) + J f{s)ds^ almost everywhere on I. 

In order to do so, the first step is to prove that ([5]) is verified in the sense of distributions. 
From <^2il it follows that 

N-l N-l 

(16) yip € V{I) , {muh,v) = rnhr<p{e) + ^ /iA"+V(i")- 

n=l n=l 

We are going to pass to the limit in this equation. By (fTSll . {miih, (f) converges to (mti, if). To 
study the first term of the right-hand side, we write 

f{s)[ip{n-if{s)]ds- / f{sMs)ds. 



hj2fMn = [ f{sMs)ds+Yl 

n=l n=l 



to 



The convergence to zero of the sum over n comes from uniform continuity of (p. Combining 
this with \t^ — t^\ = h gives 

N-l T 

V mhpp{e) — ^ m / f{s)ip{s)ds when /i ^ 0. 
n=i -^0 

The argument for the last term is similar. We write 

N-l N-l ^ti 

Xhis) [ipif") - p{s)] ds - / Xh{s)ip{s)ds. 



^ hX^+^pin = [ Xhis)p{s)ds + ^ / 

n=l n=l -^i" 



tO 
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The convergence to zero of the sum over n comes from uniform continuity of (f and lemma [2. 121 
and the last term is equal to zero for all h. This, together with lemma [2.121 gives 

Af-l 

/ir+V(*'') — > (A, f) = -(7, when h ^ 0. 

n=l 

Finally, passing to the limit in ()16p we obtain 

{mq - 7, (^) = {mf, ip) , \/ip £ V{I), 

as required. 

Then, by density of in Cq(/) and the fact that mq — 7 is in BV{I), we get 

mq — 7 = mf in M.{I). 
Integrating this equality over [0,t[ (Stieltjes integral of BV functions) we obtain 

{mq - -f){t^) - {mq - j){0~) = mf, 

Jo 

and the result follows from this, by using 7(0^) = and a.e. continuity of mq — 7. 
(4) Proof of 0^7 = 

To prove that ((7,7) is solution to (jlip . it remains to show that q^y = almost everywhere. 
For all n we have q'^-f^ = 0. However, qh'Jh is not identically equal to zero. We build new 
functions ijh and 7/1, piecewise constant, with respective values g" and 7" on ji", We 
now have qhjh = and simple computations give 

\\Qh - < Mh - Qh\\L^(I) + \\Qh - Q\\l°°(I) < h\Wh\\L°°{I) + IWh - Q\\l°°{I) 

and 

Wlh - iWl^i) < Wlh - 7h||Li(/) + Wlh - iWl^i) < 2^\^^\\Li{I) + \hh - 7||li{/)- 

Combining the first inequality with lemma 12.101 and (jl4p gives uniform convergence of ijh to 
q. Putting together the second inequality, lemma [2.121 and (fT5|) . we see that converges to 
7 in L^{I) which implies that the sequence converges up to a subsequence almost everywhere 
on /. Finally, letting h go to zero in qhjh = gives (77 = almost everywhere as required. 

This completes the proof of theorem 12.81 □ 



2.4. Validation: coupling with a fluid/particle solver. We consider the same experiment that 
the one considered in section [21 without inertia. The radius of the particle is taken equal to 1 and the 
viscosity of the fluid is fi = 3. The balance of forces reads 

(17) Vt, FiMt)) + f{t) = 0, 

where f{t) = -2 until time 2 and f{t) = 2 if t > 2. 

To obtain a reference solution we first compute, as accurately as possible, the map q Fiub,UQ (q) for a 
given velocity uq = —1 and q S [0, 1]. To do so, we begin with computing Fiub,uo{Qk) where {qk)k=i..M 
is a regular subdivision of interval [0, 1]. This is done, for each q^, solving the Stokes problem in the 
fluid with Dirichlet boundary conditions and computing the force Fiub^uoiQk) exerted by the fluid on 
the particle. The computations are carried out in tree-dimensions using an axisymmetric formulation 
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and the Finite-Element solver FreeFeni++. On the left side of figure [5l we plot the numerical results 
obtained (circles). They agree with the asymptotic expansion ([1]) for small distances (solid line). 
Finally, the map q —>■ Fiy^i,{q) is approximated using a least square approximation of the numerical 
results by a polynomial of degree 3 (dashed line). 

The reference solution is obtained discretizing the time interval and computing the velocity of the 
particle at each time-step. We write that n" = a'^UQ and, using the linearity of the lubrication 
force with respect to the velocity, we compute a" as the solution to 

«"Fz«6,uo(9") + /(n = 0. 
The trajectory obtained is plotted against time on the right side of figure [5l 




io9(q) 



Figure 5. Approximation of the lubrication force (left) and reference solution (right). 



We now want to observe the influence of the method employed to deal with contacts in fluid/particle 
simulations. To do so, we use an axisymmetrical version of the fluid/particle solver implemented 
with FreeFem++ and described in [T3]. On figure El we plot the solution given by this solver for a 
mesh size 5x = r/10 (dashed line). We can observe that the particle remains glued. Indeed, due 
to the space discretization, the characteristic function representing the rigid particle ends up with 
touching the boundary of the domain and the Dirichlet boundary condition prevents it from taking 
off. Consequently, it is necessary to deal with the problem of contact and to prevent the characterictic 
function from intersecting the boundary of the domain. Two methods are tested: the fluid/particle 
solver is coupled with an inelastic contact algorithm and with the gluey contact model. The coupling 
is performed using the splitting strategy described in remark [2.71 In each case, the constraint for the 
distance is set to q > r] with r] = 5x. The numerical results are compared on figure El We observe 
that, for the inelastic model (solid line with crosses), the particle takes off as soon as it is pulled. To 
the contrary, using the gluey contact model (solid line with circles), the particle remains glued and 
the trajectory finally joins up with the reference one. This is a validation of the gluey particle model 
and it emphasizes the necessity to take the lubrication force into account when dealing with contacts. 
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t 



Figure 6. Comparison of the numerical solutions for different contact models. 

Finally, we observe on figure [7| the behaviour of the two contact models with respect to the parameter 
r] which is the minimal distance allowed between the particle and the plane. We can see that the 
trajectories obtained for different r] separates after unsticking time when using the inelastic contact 
model (left side of the figure). This is due to the fact that, for this model, the particle unsticks as 
soon as it is pulled. To the contrary, the gluey particle model is not so sensible to parameter r/ (right 
side of the figure). 




Figure 7. Impact of rj on the numerical solution for inelastic (left) and gluey (right) contact. 



2.5. Extension to rough solid surfaces. As already said, it has been proved that smooth solids 
can not undergo contact. However, from our experience, we know that the particle should touch the 
plane in finite time. One of the reasons explaining this behaviour is that physical particles are not 
smooth. Recent experiments described in [261 show that the lubrication force exerted on a rough 
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particle is the one would be exerted on a shifted smooth particle: 




Figure 8. Equivalent smooth sphere. 

contact with the plane {q + Qs doesn't go to zero in finite time) but the real surfaces can collide {q can 
go to zero). 

Taking these results into account in the gluey particle model, we consider that, as soon as g = ri^s+'^2,s 
(see notations on figure[9]), there exists a real solid/solid contact. During this contact, the forces acting 
on the particle are not registered anymore. To model such a behaviour, it suffices to recall that 7 is 
the limit of 7^ = 67r/iZn(g^) (see remark [2. 4p and to impose 

7 > 67r/x/n(rs_i +rs,2)- 

The trajectory computed for this model is plot on figure [H We can observe that the rough particle 
takes off before the smooth one. Contrary to what has been said in remark 12.41 for the smooth case. 




Figure 9. Rough solids: notations (left) and gluey particle model (right), 
it is now important to know the value of 7 in order to truncate it. Therefore, it is essential to take 
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the radius into account in its evolution and to use equation Q in the gluey particle model: 



In that case, the trajectory of the particle depends on r. 

From an algorithmic point of view, this rough model can easily be taken into account in algorithm 12.61 
by changing step ([3|) in 

7"+i =7"- Aa"+i, 
where r is the radius of the sphere and by adding the following step: 

3. Multi-particle case 

3.1. Modelling. We generalize the gluey particle model Q-® to the multi-particle case. We consider 
a system of N spherical particles in three-dimensions. Xj stands for the position of the center of particle 
i in IR^ and fj G IR^ for the external force exerted on it. Let x G IR^^ be defined by x = (..., Xj, .. .) and 
f € IR'^^ by f = (. . . , fj, . . .). We denote by Dij the signed distance between particles i and j, and Sij 
by eij(x) = (xj — Xj)/||xj — Xj|| (see Fig. llOp . We define M as the mass matrix of dimension 3A^ x 3A^, 
M = diag{. . . ,mi,mi,mi, . . .). Vector Gij € IR^''^ is the gradient of distance Dij with respect to the 
positions of the particles: 

Gij(x) = VxAi(x) = (...,0, -ei_,(x) ,0,...,0, eij(x) ,0,...,0)*. 

i 3 




Aj(x) 



Figure 10. Particles i et j : notations. 



In that context, there are N{N — l)/2 pair of particles and we denote by 7 = (. . . ,7jj, . . .) G |R-^(^ i)/^ 
the associated sticking variables: "jij is stricly negative if particles i and j are glued. Then, using the 



fact that ^ ^•'^ 



dt 



Gi,(x) 



X, we define the following space of admissible velocities: 



Gij(x)-V = Oif -fij{t-)<0 

G^j{x) • V > if 7ij(r) = 0, Dij{t) = 



To finish with notations, we denote by A = (. . . , Xij, . . .) G |R^(^ the vector made of the Lagrange 
multipliers associated to these A^(A^ — l)/2 constraints. 
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The multi-particle model is the natural counterpart of the particle/plane one: 

Mx = Mf + ^ AijGij(x), 

(19) < supp(Aij) C {t? Aj(t) = 0} for ah 
7 = -A, 

Aj > , 7ij < for ah i, j, 

_ x(0) = x° St. Aj(0) > for ah i,j , i(0) = u° , 7(0) = 0^n(n-^i)/2. 

Remark 3.1. The Lagrange multiplier Ay, associated to the constraint between particles i and j, is 
activated (non zero) only if these particles are in contact. The additional force due to this contact is 
equal to AjjGij(x). From the expression of Gij(x), we get that this force only concerns the particles 
involved in the contact: it is equal to — Aijejj(x) on particle i and Ajjejj(x) on particle j. 

Remark 3.2 (Roughness and radius). As for the particle/plane case (see section [23]) . roughness can 
be taken into account by imposing a threshold on 7: 

6Trfiln{ri^s + rj^s) < lij for all i, j, 

where n is the size of roughness of particle /. As noticed in the particle/plane case, it is now important 
to take the radius of the particles into account in the evolution of 7. To do so, in the same way as 
in the particle/plane case, we come back to the way the gluey particle model has been built and take 
into account all the constants involved in the first order assymptotic developpement of the lubrication 
force exerted between two particles We obtain the following evolution equation for 7: 

7 = -R\ 

where R is the diagonal matrix of dimension N{N — l)/2 with coefficients Rij,ij = {vi + 'fj)'^/{rfr'j). 



3.2. Algorithm. Let h be the time step. We denote by V" = (..., V", . . .) G IR^^ the approximated 
velocities of the particles at time = nh. Let x", 7" and A" be the respective approximations of x, 
7 and A at time t". 



The discretization of the continuous constraints Cx,-y(t") is inspired by [18] and corresponds to a first 
order approximation of the constraints: 



K(x",7") 



V G R^^ s.t. 



Aj(x'^) + /iGi, (x") • V > if 7-. = 
A,(x") + /iGi,(x") • V = if 7". < 



Using this discrete space of admissible velocities, the time discretization of (jl9p is now a direct adap- 
tation of algorithm 12.61 to the multi-particle case. 
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Algorithm 3.3 (Multi-particle). For all n > 0, let x", V^, 7" and A" be given. We define P 



iTl+l 



f(s)ds. 



(1) Computation of the a priori velocity, without taking the lubrication force into account 

yn+l/Z ^ Y" + /if". 

(2) Projection of the a priori velocity on the set of admissible velocities, 



Y'"+i _ ■yn+1/2 



= mm — 

M ve-ft'(x",7") 2 



V - v"+^/2 



M 



From this projection step, we obtain A"^"*^. 
(3) Updating of-f. 



n+l 



0. 



(4) Updating of x, 



X 



n+l 



x" + hV 



n+l 



Remark 3.4. In the same way as in section 12.51 and remark 12.71 for the particle/plane case, this 
algorithm can be extended to rough solids and coupled with fluid/particle solvers using a splitting 
strategy. 

Remark 3.5 (Obstacles). Suppose there exists A'^o fixed obstacles (walls of a box containing the 
particles for example). It is straightforward to add the NNq new constraints in i^(x",7"). Now, 
suppose these obstacles are moving with a prescribed velocity. We denote by y""*"^ the (known) vector 
giving their position at time t""*"^. The space of admisible velocities becomes: 

Pairs (i,j) particle/particle : 

Ai(x") + /iGij(x") • V > if 7f^. = 

Ai(x") + hGijix.'') • V = if 7;;. < 

Pairs {i,k) particle/obstacle : 

Afc(x", y"+i) + /iGifc(x", y"+i) • V > if 7,"^ 



/^(x-,y"+\7"^ 



V E s.t. 



> . 







Afc(x",y"+^) + /iGifc(x-,y-+^) • V = if 7^^ < J 



3.3. Finding neighbours. The most time consuming step in algorithm l3.3l is the projection step 
It is performed using a Uzawa algorithm which imposes to run matrix/vector products involving the 
contacts. However, in order to simulate large collections of particles, it is essential to avoid loops over 
the N{N — l)/2 possible contacts. To do so, we notice that it is not necessary to take into account 
all contacts at each time-step. Indeed, two particles i and j far enough to each other at time t" won't 
stick at time i""*"^ and consequently, the corresponding constraint won 't be activated (ie. A"^+^ = 0). 
We denote by Dneigh the distance above which we consider that two particles are not likely to touch 
next time-step. Then, the set of pairs of particles one has to consider at time i" is: 

C„eig/i(x") = {{i,j) G [l,Nf , i<j and Ai(x") < -D„eigh} • 
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If the pair is in the set C^ieighi^'^) ^ ^6 say that particles i and j are neighbours. Two particles 
that are not neighbours at time t" won't stick at time t""*"^ and consequently, one can restrict the set 
of constraints at time i" to: 

V(i,i) e C7„eig/i(x") , 



V G IR^^ s.t. 



A,(x") + /iG,,-(x") 
A,(x") + /iGy(x") 



V > if 7^- > 

V = if 7,";- < 



Remark 3.6. This idea not to take into account particles far away from each other is generally 
used when considering particles interacting through near field interaction forces, decreasing with the 
distance. In that case, it consists in considering that the force is negligible above a certain distance 
and consequently, it is an approximation of the model. In our case, no approximation is made. Indeed, 
if Dneigh has been chosen sufficiently large, we know that the pairs of particles that are not belonging 
to C„eig/i(x") won't interact at time t'"'^^. For example, we can choose a time step in order to limit the 
displacement of the particles to twice their radius and then set the value of Dneigh to a few radiuses. 



To construct Cneighi^"') avoiding the computation of the A^(A^ — l)/2 distances, we choose a bucket 
sorting type algorithm. It consists in dividing the computational domain into boxes of size > Dneigh 
and to compute distances only for pairs of particles belonging to neighbouring boxes (see Fig. [TT]) . 
Note that, because of step (l3|), it is not sufficient to erase at each time-step the former set of neighbours 
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Figure 11. Algorithm to find neighbours: neighbouring boxes and distances actually computed. 

and to create the new one: one has to transfer the value of 7[J- if particles i and j are in contact during 
these two successive time steps. 



3.4. Object Oriented Programming Method. To build this code, we chose to use the object 
oriented programming method for mathematical problems CsiMoon [12]. As a consequence, both 
numerical methods and models can be easily changed. For example, new methods can be chosen and 
added to the code in order to perform the projection step and to construct the set of neighbours. This 
programming method also alows us to take into account various models of external environment (dry 
environment, fluid, obstacles of different shapes...), of interparticular interactions (cohesion force...) 
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and of contacts (inelastic, gluey model, aggregation...). This leads to a modular C++ code SCoPI [T5] . 
allowing Simulations of Collections of Interacting Particles. This code has already been used to 
simulate gluey particles, crowd motion, wet particles and red-cells (as an assembly of rigid particles). 

4. Numerical simulations 

We present in this section numerical simulations of collections of gluey particles. For visualization 
reasons, we only propose here two-dimensional simulations: even though the code is intrinsically 
three-dimensional, the motion of the particles is restricted to a vertical plane. These simulations 
demonstrate that the algorithm enables to take great numbers of gluey particles into account. This, 
together with section \TM shows that coupling the gluey particle algorithm with fluid/particle solvers 
will make it possible to simulate dense fluid/particle flows, taking the lubrication force into account 
with accuracy. 

4.1. Gluey lotto: influence of roughness. The aim of this simulation is to observe the influence 
of roughness on the behaviour of multi-particle systems governed by the gluey particle model. We 
consider a two-dimensional "gluey lotto" made of 160 particles in a squared rotating mixer operator. 
The side lenght of the box is 0.5 and the radiuses of the particles are taken between 0.007 and 0.015. 
All particles have the same mass m = 1 and the gravity constant g is taken equal to 10. The 80 
particles initially situated in the left compartment of the box are black and the 80 other ones are 
white. We represent side by side on figure [12] the configurations obtained at different time steps 
for 7mm. = on the left (inelastic contacts), 7mm = — 1 in the middle (gluey rough particles) and 
7mm = —CO on the right (gluey smooth particles). In case of smooth particles, the heaps of particles 
take off from the wall when they are at the top of the box: as suggested by the particle/plane model, 
they take off only when the gravity has balanced the forces it has itself exerted to push the particles 
on the bottom wall. In the rough case, they take off earlier. 

4.2. Sedimentation of 3000 gluey particles. We consider 3000 gluey particles sedimenting under 
gravity with radiuses between 0.015 and 0.025. They are initially situated above a funnel (random 
sample of positions) with velocity equal to zero. All the particles have the same mass m = 2 and the 
gravity g is taken equal to 10. Below the funnel, a wheel rotates around its axis with angular velocity 
uo = —2 and throws the particles on a leaning fixed plane situated below it. Then, the particles slip 
along the plane and finally fall in a container. Some spherical obstacles of radius r = 0.1 are fixed on 
the plane to slow the particles movement. A threshold is imposed on 7 (7 > —10) to model roughness. 
Snapshots of this simulation are presented for different time-steps on figure [131 The code also allows 
us to model dry granular flow involving inelastic contacts. In figure [H] we compare the configurations 
obtained at the same time-step for such a simulation and the previous gluey one. Finally, we plot on 
figure [l5] the values of 7 for a given configuration of the gluey simulation. For each contact, a tube is 
plotted between the two involved particles and, the larger is 7,^ (ie. the more the particles are glued), 
the more the grey is dark. We can see the network of the forces leading to a packed configuration in 
the funnel. The particles are smoothly unsticking from each other when leaving the wheel. 
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Figure 12. Gluey lotto: configurations at different time-steps for 7mm = on the left 
(inelastic contacts), 7mm = — 1 in the middle (gluey rough particles) and 7mm = — oo 
on the right (gluey smooth particles). 




Figure 13. Snapshots of a two-dimensional simulation, 3000 particles: configurations 
at time-steps n = 8826 - 14226 (top) and n = 18726 - 30726 (bottom). 
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Figure 14. Snapshot of two-dimensional simulations, 3000 particles: dry (left) and 
viscous (right) simulations at time-step n = 14226. 
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Figure 15. Snapshot of a two-dimensional gluey simulation, 3000 particles: configu- 
ration and values of 7 at time-step n = 14226. 
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